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Full seven-equation Reynolds stress turbulence models are promising tools for today’s 
aerospace technology challenges. This paper examines two such models for computing 
challenging turbulent flows including shock-wave boundary layer interactions, separation and 
mixing layers. The Wilcox and the SSG/LRR full second-moment Reynolds stress models have 
been implemented into the FUN3D unstructured Navier-Stokes code and were evaluated for 
four problems: a transonic two-dimensional diffuser, a supersonic axisymmetric compression 
corner, a compressible planar shear layer, and a subsonic axisymmetric jet. Simulation results 
are compared with experimental data and results computed using the more commonly used 
Spalart-Allmaras (SA) one-equation and the Menter Shear Stress Transport (SST-V) two- 
equation turbulence models. 


Nomenclature 


a = Reynolds-stress anisotropy tensor 

b = mixing layer thickness 

skin friction coefficient 

coefficients of the Speziale-Sarkar-Gatski model 
coefficients of the Wilcox RSM pressure-strain 
= equilibrium parameter 

Reynolds stress diffusion coefficient 

= specific diffusion tensor 

= distance to nearest wall point 

Mentor’s blending function 

local diffuser height 

= diffuser throat height 

specific turbulent kinetic energy 

turbulent mass flux tensor 
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= specific turbulence production tensor 
= static pressure 

radial coordinate* 

= strain rate tensor 
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traceless strain rate tensor 

= time 

= velocity component in the x-direction 
local mean streamwise velocity 
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freestream velocity difference, AU = U, — U, 
velocity component 

rotation tensor 

= spatial coordinate in the streamwise or axial direction 
= Cartesian coordinates 

= spatial coordinate in the vertical or wall-normal direction 
= coefficient of @ production 

= closure coefficient 

= coefficient of @ production 

coefficient of m destruction 

= closure coefficient 

= closure coefficient 

= coefficient of destruction 

closure coefficient 

= Kronecker delta 

turbulent dissipation rate 

= argument to F, blending function 

= dynamic viscosity 

= specific pressure-strain correlation tensor 

= density 

= coefficient of w diffusion 

= closure coefficient 

coefficient of cross-diffusion 

Ow = coefficient of w diffusion 

Ti = viscous stress tensor 
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¢© = coefficient of ¢ equation 
¢™ > = coefficient of w equation 
w = specific dissipation rate 


= fluctuating component, d"” = ¢—-—¢ 
= Reynolds averaged component, o = Jim Pa peotre pdt 
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T = turbulent 

t = total 

0 = stagnation condition 
co = freestream condition 


I. Introduction 


OR aerospace propulsion flows, the most common practice in computational fluid dynamics (CFD) analyses is to 

use Reynolds-Averaged Navier-Stokes (RANS) codes with one- or two-equation turbulence models. Current 
RANS turbulence models predict steady, fully turbulent attached flows at all speed regimes reasonably well, but are 
still unable to reliably predict most separated flows. Large-eddy simulation (LES) and direct-numerical simulation 
(DNS) methods are being used for some applications, however they require very fine grids for wall bounded flows 
and at high Reynolds numbers, and therefore will not be practical for many years. [1, 2] Hybrid unsteady RANS/LES 
methods are increasingly common for certain classes of simulations including separated flows, although the techniques 
to combine the near-wall RANS region with the outer, large-eddy simulation region need further development. [1] 

For aerospace propulsions flows, RANS will be used for a significant portion of CFD analyses for the foreseeable 
future, due to limitations in computational power. [1] Traditional RANS linear and nonlinear one- and two-equation 
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turbulence models, such as Spalart-Allmaras (SA) [3] and Menter Shear-Stress Transport (SST) [4], are frequently 
used. These linear models incorporate the Boussinesq approximation to give the Reynolds shear stress tensor in terms 
of the mean strain rate tensor and the eddy viscosity. [5] Nonlinear models incorporate additional higher-order terms 
that are functions of the mean strain and rotation rate tensors. An alternative to these approaches is to use a more 
advanced form of RANS turbulence modeling, where the individual transport equations for the Reynolds stresses are 
solved, allowing for a more detailed representation of the flow physics. These models are known as full second- 
moment Reynolds Stress Models (RSMs). 

This paper describes a study using the FUN3D [6] unstructured CFD code to evaluate two RSMs: the Wilcox 
RSM [5] and the SSG/LRR RSM [7]. This work supports the NASA Revolutionary Computational Aerosciences 
technical challenge to, “Identify and down-select critical turbulence, transition, and numerical method technologies 
for 40% reduction in predictive error against standard test cases for turbulent separated flows, evolution of free shear 
flows and shock-boundary layer interactions on state-of-the-art high performance computing hardware.” [1] The 
first step in this study was to reproduce previous results obtained with the Wind-US code using the SA and SST 
turbulence models in order to confirm that the turbulence models are implemented equivalently in each code. Next, 
FUNS3D was run with the two RSMs. Four cases were examined: a transonic diffuser [8, 9, 10], supersonic flow over 
an axisymmetric compression corner [11, 12, 13], a compressible planar shear layer [14], and a subsonic axisymmetric 
jet, [11, 15]. Results were compared with the SA and SST models, and with experimental data. 


II. The Codes 


The FUN3D code was used for the computations described herein and compared with results previously computed 
using the Wind-US code. Both are NASA production codes using Reynolds-averaged flow solvers and are described 
briefly below. 


A. Wind-US 

The Wind-US code [16, 17] has both structured and unstructured solvers. The flow equations are evaluated using 
second-order-accurate finite differences defaulting to second-order Roe for structured grids and second order HLLE 
for unstructured grids, although other schemes may be specified in the user inputs. The partial differential equations 
are written in conservative form and explicit terms are modeled using either upwind or central differencing. Implicit 
terms are computed using either an approximately factored or four-stage Runge-Kutta scheme. 

The Wind-US results described in this paper are previously reported computations using the well-known Spalart- 
Allmaras (SA) [3] one-equation and Menter shear-stress transport with vorticity source term (SST-V) [4] two-equation 
turbulence models. The results are compared with FUN3D results computed using the same turbulence models in 
order to verify that these baseline turbulence models are implemented equivalently in each code. Additional details 
can be found at www.grc.nasa.gov/www/winddocs and in Ref. [18]. 


B. FUN3D 

FUN3D [6] is a node-centered, unstructured implicit solver developed by researchers at the NASA Langley 
Research Center. It uses finite-volume discretization and is formally second-order accurate in space. Explicit terms 
are calculated using Roe’s flux difference splitting, however other methods are available. More information about 
FUN3D can be found at fun3d.larc.nasa.gov and in Ref. [6]. Note that at the time of this study, the full Reynolds 
stress models in FUN3D cannot be used for periodic grids whose side planes are rotated through a small angle. Instead, 
it requires 90-degree grids with both axes aligned with one of the constant x-, y- and z-coordinate surfaces. For the SA 
and SST-V models, small-angle periodic boundary conditions are available, so grids one-cell wide in the 
circumferential direction were used. 


C. Solver Inputs 

The inputs to both codes were set at the values suggested for each case according to the code documentation or 
recommended practices. For Wind-US, the right-hand-side inviscid explicit operator was set to the second-order 
HLLE flux splitting algorithm for the transonic diffuser case which was computed using an unstructured grid and the 
second-order Roe upwind scheme for the other cases, which were computed using structured grids. For FUN3D, the 
Roe second order upwind scheme was used. All cases were run steady state, with the exception of the axisymmetric 
jet, which required a time-accurate solution to reach convergence. 
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Ill. The Turbulence Models 


The full RSM turbulence models evaluated in this work are the SSG/LRR RSM and the Wilcox RSM. The more 
standard Spalart-Allmaras (SA) one-equation model [3] and the Menter shear stress transport (SST-V) two-equation 
model with a vorticity source term [4] were also used for comparison with the RSMs. Details about all of these models 
are available on the Turbulence Modeling Resource (TMR) website, turbmodels.larc.nasa.gov [19]. For one case, the 
mixing layer flow, a nonlinear explicit algebraic stress model (EASM), available in Wind-US, was also used. [20] 
The SA and SST-V models use the Boussinesq approximation, 


. « ls 2 (1) 
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a constitutive relation, to compute the turbulent stress tensor, T; pp defined below as, 
Tij “ ee) (2) 


On the other hand, the RSMs solve transport equations for each of the six unique Reynolds stresses. A seventh equation 
is required to determine the length scale variable, w. 

The SSG/LRR full Reynolds stress model, as described on the TMR, and in Refs. [7] and [21], was developed 
under the European Union project FLOMANIA and consists of the six equations of the Reynolds stress transport 
equations, plus Menter’s baseline @ equation for the length scale. The Reynolds-stress transport equation is given by 


Eq. (3). 
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The blended Speziale-Sarkar-Gatski/Launder-Reece-Rodi pressure strain model as defined in Refs. [7, 19] is given 
by Eq (4). 
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The remaining terms in Eq. (3) are given below. 
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The turbulent mass flux term, M;;, is assumed to be negligible. 


jy? 


A length scale equation is required to close the system. The SSG/LRR model uses an w-equation similar to 
Menter’s SST-V model which blends an w-equation near the wall with an ¢ equation in the outer boundary layer. 
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The isotropic dissipation rate is defined as 


e=C,ko (13) 


with C, = 0.09. The equations used for blending the coefficients 6 = a, B,,0,0q are given below. 
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where d is the distance to the nearest wall. The pressure-strain coefficients are blended between Launder-Reece-Rodi 
(LRR) [22] near walls (without wall-correction terms) and Speziale-Sarkar-Gatski (SSG) [23] away from walls. The 
coefficients are given in Table 1. 


Table 1. Coefficients for the SSG/LRR full Reynolds stress turbulence model. 
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The diffusion term given in Eq. (11), above, is a generalized gradient diffusion model [24]. In the FUN3D input 
file, the SSG/LRR RSM is specified as “SSG/LRR-RSM-w2012.” There is also a version of the SSG/LRR RSM 
which uses a simple diffusion model [25] and is specified as “SSG/LRR-RSM-w2012-SD” in the input file. In this 
model, the diffusion term is modeled as given in Eq. (17) 


: a[/_ D_ \a(t) eee 17 
pDi; = — ax, (« — C, ur) ax, Simple Diffusion ay) 


with: 


2 
D=0.5¢,F; + 30.22(1 — F,) (18) 


Both diffusion terms are used in the paper. The transonic diffuser and mixing layer simulations both used the 
standard diffusion model. The simple diffusion model, which tends to be more stable, was used for the shock-wave 
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boundary-layer interaction and the axisymmetric jet, because these simulations had more difficulty reaching 
convergence. 

The Wilcox RSM as described on the TMR and in Ref. [5] also solves the six equations for the Reynolds stress 
tensor and an omega equation for the specific dissipation rate, w. The Reynolds stress equation is given by Eq. (19). 
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The specific dissipation rate is given by Eq. (20). 
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The pressure strain correlation is given by: 
a 2 fe 2 . 2 2 1 


with, 


1 a 
P= 5 Phx and Ur = pk/w (22) 


The production, P;;, given by Eq. (5), the simple strain rate tensor, S;;, given by Eq. (8), and the averaged 


rotation tensor, W;;, given by Eq. (10), and are used in both the SSG/LRR and Wilcox RSMs. The closure 
coefficients for Eq. (21) are given in Table 2, and in Eqs. Eqs. (23) and (24), below. 


Table 2. Closure coefficients for the Wilcox full Reynolds-stress turbulence model. 
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The EASM [20], available in Wind-US, was used for comparison with the other models for the mixing layer case. 
This model is derived from a reduced form of the Reynolds Stress transport equations. The resultant expression for 
the Reynolds stresses is similar to Eq. (1), but includes terms that are non-linear in the strain and rotation rate tensors. 
This particular EASM uses a simplified form of the SSG pressure-strain relation that is a linear function of the 
anisotropy. EASMs provide more information about the shear stress than the SA and SST-V models, while remaining 
less computationally expensive than the RSMs. Wind-US EASM results are shown for the mixing layer case in section 
IV(C). 
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IV. Test Cases and Results 


A. Transonic Diffuser 

1. Experimental Setup 

The first case examined was transonic flow through a converging-diverging diffuser. This flow has a strong normal 
shock wave in the throat causing a separated region to form on the top wall. This case was run previously with Wind- 
US, and was used to verify that the SA and SST-V turbulence models produce the same result in both Wind-US and 
FUN3D. Once this verification was complete, the flow was computed using the RSMs. Comparisons were made with 
the experimental data obtained by Sajben et al. [8] The duct geometry in the x-y plane is shown along with Mach 
contours computed using the SA model in Figure 1. The geometry is a convergent/divergent channel with a flat 
bottom and a contoured top wall. The flow conditions at the inflow plane are Mach 0.9 with a total pressure and 
temperature of 19.58 psi and 540 degrees R. A normal shock wave forms downstream of the throat at x/Hrproat = 
2.24, causing the flow on the top wall to separate. It re-attaches at approximately 6 throat heights downstream. 

2. CFD Methodology 

The unstructured grid used for this case, shown in Figure 2, contains hexahedral cells near the walls and tetrahedral 
cells in the center. A constant area section extends 10 throat-heights downstream of the diffuser exit station to 
eliminate any boundary layer separation at the outflow as the flow sets up. (Only a portion of this is shown in Figure 
2.) The grid contains 54,854 points, with 301 points in the axial direction. In the vertical direction, there are 41 
vertical grid points in the hexahedral grid near each wall and approximately 91 total points in the vertical direction. 
The grid is one throat height deep in the transverse direction. The y* values one point from the wall are approximately 
1.0. This grid is also described on the NPARC Alliance Verification and Validation Website [26]. Results computed 
on this grid were equivalent to those computed on a grid twice as dense. 

The boundary conditions used are those described in Ref. [26] to set up a strong shock in the diffuser. The top 
and bottom walls are adiabatic and viscous. The inflow conditions are set to the Mach 0.9 conditions described above. 
The outflow static pressure is set to 14.775 psi to allow the normal shock to form at the desired location in the throat, 
X/Henroat = 1.98. The full diffusion version of the SSG/LRR RSM was used. 


4 3 -2 0 1 2 3 “ 5 6 7 8 


X/Hrhroat 


Figure 1. Transonic diffuser. Mach contours computed with Wind-US and the SA turbulence model. 
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Figure 2. Transonic diffuser. Computational mesh. 

3. Results 


The results computed using Wind-US and FUN3D with the SA and SST-V models are shown in Figure 3 and 
Figure 4. The normalized pressure on the top and bottom walls in Figure 3 indicates that both codes give equivalent 
results. Both codes and turbulence models predict the pressure rise at the normal shock wave. However, just 
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downstream of the shock, at x/Hnyoa¢ between 3 and 6, all of the results slightly overshoot the experimental values. 
Further downstream, they agree with the experiment. The velocity profiles shown in Figure 4 also indicate equivalent 
results between the two codes. At x/Htnroat = 2.9 and 4.6, the CFD results are all correctly predicting that the flow 
is separated on the top wall, and the velocity profiles are matching up well with the experiment. At x/Htnroat = 6.4, 
however, the CFD results still indicate that the flow on the top wall is separated whereas the experiment indicates that 
it is attached. At x/Htnroat = 7.5, the flow is attached for both the CFD and experiment, however, the CFD results 
predict thinner boundary layers on the top and bottom walls. 
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Figure 3. Transonic diffuser. Normalized wall pressure on the bottom and top walls computed with Wind-US 


and FUN3D. 
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Figure 4. Transonic diffuser. Velocity profiles computed with Wind-US and FUN3D. 


The FUN3D results computed with the Wilcox RSM and the SSG/LRR RSM are plotted with the SA and SST-V 
results in Figure 5-Figure 7. In the pressure plots of Figure 5, the RSM results are slightly better at predicting the 
pressure rise just downstream of the shock, however, the values further downstream are slightly higher than the 
experimental values, whereas the SA and SST-V results agree with the experimental values for x/Hinyogt values 
between 6.5-9. The velocity profiles are shown in Figure 6. The RSMs do not predict separation, but they do a better 
job of predicting the boundary layer shape at x/Htnroat = 6.4 and 7.5. The RMS u-velocity is shown in Figure 7. It 
was computed by taking the square root of the uw’ from the FUN3D output. For the SA model, w'u' is computed 
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from the Boussinesq approximation (Equation (1)) using Bradshaw’s approximation for k, where k = 
Pref S55, /(0.31p). The plots show that the RSMs are better than the SA and SST-V models at predicting the 
values at the peaks of the profiles. For x/Htnroat = 2.9 and 4.6, the values at the center of the diffuser are under- 
predicted for all of the models. 
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Figure 5. Transonic diffuser. Normalized wall pressure on the bottom and top walls computed with FUN3D. 
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Figure 6. Transonic diffuser. Velocity profiles computed with FUN3D. 
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Figure 7. Transonic diffuser. RMS u-velocity computed with FUN3D. 


Overall, it is unclear if the RSMs offer benefits over the traditional turbulence models based solely on this case. 
For the pressure, they are better at predicting the pressure rise just downstream of the shock, but then over-predict it 
further downstream. For the velocity profiles, the RSMs are worse than the traditional models at predicting the 
separation on the top wall, but give better agreement with data further downstream. For the 
Urms profiles, the RSMs are significantly better at predicting the profiles. 


B. Axisymmetric Supersonic Compression Corner 

1, Experimental Setup 

This study examines supersonic flow over an axisymmetric compression corner [11, 13, 15, 27, 28]. The model 
is a 5.08 cm diameter cylinder with a 30-degree flare, which generates a shock wave (Figure 8). The cylinder has an 
upstream cusped nose designed to minimize the strength of the shocks, and data indicates that reflected shocks from 
the tunnel walls have no effect in the measurement region of interest. The flare is located 1.0 m downstream of the 
cusp-tip, allowing a turbulent boundary layer to develop upstream of the shock-wave boundary-layer interaction. The 
flare surface begins at x = 0.0 cm and ends at x = 5.196 cm. The test section has a Mach number of 2.85, a Reynolds 
number of 16x10°%/m, a stagnation pressure of 1.7 atm and a stagnation temperature of 270 K. The data consists of 
laser Doppler velocimeter mean velocities and Reynolds stresses, surface static pressures, Schlieren photography, oil 
flow visualizations, and holographic interferometry data. 

This flow has a complex shock-wave/boundary-layer interaction region, as shown by the FUN3D results in Figure 
9, depicting Mach contours in the interaction region computed with the SA model. A primary shock wave is generated 
when the flow encounters the flare surface. The shock penetrates the subsonic portion of the boundary layer causing 
a strong adverse pressure-gradient, which leads to separation in the corner between the cylinder and the flare surfaces. 
The separation extends upstream of the flare, inducing a secondary shock. On the flare surface, the flow at the 
reattachment point is turned parallel to the flare, generating a shock that coalesces with the primary shock. The flow 
is turned back axially downstream at the cone-afterbody, causing the flow to expand and weaken the shock wave. 
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Figure 9. Axysimmetric compression corner. Mach contours computed using the SA turbulence model. 


2. CFD Methodology 

For this case, a single-cell, axisymmetric wedge grid was used for computations with the SA and SST-V models. 
Since wedge boundary conditions are not yet available in FUN3D for the RSMs, a 90-degree version of the grid was 
generated for use with the RSMs. The axisymmetric grid was originally generated by Debonis [12] using Pointwise 
[29] for structured Wind-US computations. It begins on the cylinder body, downstream of the cusp. The length of 
the cylinder, 75 cm, was chosen such that the CFD matched the experimental boundary-layer profile at the first 
measurement station. The grid was thoroughly tested for grid convergence in Ref. [12], resulting in a dense grid in 
order to capture the shock and separation features. The grid was orthogonal to the wall, with wall spacing of 
5.0x10° cm and a maximum y* of 0.2. Grid lines are parallel to the shock location with clustering near the expected 
shock location. There are a total of 1,265 axial points and 729 radial points for a total of 922,185 points. Pointwise 
was used to save the axisymmetric grid in unstructured, hexahedral format compatible with FUN3D. For the 
computations using the RSMs, the axisymmetric grid was extruded 90 degrees about the x-axis at 5.6 degree 
increments for a total of 17 equally-spaced points in the circumferential direction and a total of 15,478,857 points. 
The grid is shown in Fig. 10. 

The FUN3D boundary conditions were set equivalent to those used by DeBonis for Wind-US. The inflow 
boundary was fixed using the supersonic freestream conditions. The outflow was extrapolated. The upper boundary 
was a farfield condition which used Riemann invariants. No-slip adiabatic viscous conditions were set at the solid 
walls. The Roe inviscid flux method [30] with the min-mod inviscid flux limiter [31] was initially used for all of the 
cases. The SST-V model had difficulty reaching a converged solution, where solutions were considered to be 
converged with the most sensitive parameters of the solution, the skin friction and the shear stress profiles, unchanging. 
FUN3D was first run with the SA and SST-V models, and results were compared with the Wind-US results of 
DeBonis. The SA results were in good agreement between the two codes, confirming that the FUN3D inputs and 
boundary conditions were set to appropriate values. The SST-V solution did not reach convergence, so the flux limiter 
was set to Venkatakrishnan [32] instead of min-mod. This gave results with converged solution variables that were 
in good agreement with the Wind-US SST-V results. However, the residuals showed more oscillation than desired. 
The solution was then run using the Van Albada flux limiter with a heuristic pressure limiter [33] and with a 
selectively-dissipative version of Edwards low-dissipation flux splitting scheme (LDFSS) [34]. Freezing the flux 
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limiter after 1000 iterations allowed for the mean flow and turbulent residuals to reach close to machine zero levels. 
The resulting solution had smoother residuals. However, the onset of separation was further downstream and the 
pressure rise was too high. These results are not included in the plots, as they most likely reflect the effect of the 
dissipation of the numerical methods, rather than the effect of the turbulence model. The simple diffusion version of 
the SSG/LRR (SSGLRR-w2012-SD) model was used. These results also had oscillatory residuals, but the solution 
stopped changing after 30000 iterations, and so the results are included for comparison. It must be emphasized though, 
that without full iterative convergence, these results are incomplete, due to the equations of motion not being entirely 


satisfied. 
Me i 


oy 


i 
LY 
Hy g 
4 
is 
MH 
Me 
ae 


CRATER 


(b) Close-up of flare, every 8" point shown. 


Fig. 10. Axisymmetric compression corner. Grid. 


3. Results 
For this case, we examined the static pressure and skin friction coefficient on the model, axial and radial velocity 


profiles, and shear stress profiles. The surface static pressure is shown in Figure 11a. The first pressure rise upstream 
of the compression corner is caused by the shock-wave induced corner separation. The larger pressure rise beginning 
at the compression corner (x = 0 cm) is caused by the shock wave induced by the flare. The pressure increases until 
it reaches the end of the flare, then drops as the flow expands due to the flow turning back axially. The SA model 
shows the best agreement with the experimental pressure values. It predicts the start of the first pressure rise 
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accurately, and predicts the correct pressure rise at the corner and does well on the flare. The small blip in pressure 
at the corner is a result of a stagnation region caused by a counter-rotating vortex in the compression corner. The 
SST-V and the Wilcox RSM behave very similarly and predict the first pressure rise at approximately x = -4.5 cm, 
which is upstream of the experimental value. The SST-V model under-predicts the maximum pressure rise on the 
flare slightly, while the Wilcox RSM has a slightly higher overall pressure rise, similar to SA. The SSG/LRR RSM 
predicts the onset of separation much too far upstream, at approximately x = -6 cm, and under-predicts the pressure 
on the flare region. 

The skin friction coefficient, C;, is plotted in Figure 11b and indicates the attached and separated regions in the 
flow. Upstream of the separation, the reported range of experimental Cy values is 0.00155-0.0017, [28] and the 
separation begins between -3.25 and -2.75 cm [12]. All turbulence models predict the correct skin friction coefficient 
values upstream of the separation. The separation locations correspond with the pressure rise of Figure 11a. The SA 
model solution shows the separation beginning at x = -2.78 cm, which agrees with the experiment. The SST-V model 
and Wilcox RSM predict the onset of separation at similar locations upstream of the experimental value, at -4.40 cm 
and -4.35 cm, respectively. The SSG/LRR RSM predicts beginning of the separation even further upstream at -5.74 
cm. All of the models show a rise in Cy, in the corner due to the counter-rotating vortex, with the SSG/LRR model 
having the largest vortex. The Wilcox and SSG/LRR RSMs reattach at x=1.94 and 2.05 cm, respectively, and the SA 
and SST-V models reattach a little further downstream at 2.34 cm and 2.44 cm. The spike in Cy, at the aft end of the 
flare, is due to the high shear stress as the flow accelerates and expands around the corner. The separation and 
reattachment values are summarized in Table 3. 


02 6 Experiment 
SA 
01g |———— SST-v 
————— SSG/LRR RSM 
0.16 ——— Wilcox RSM 


4 6 4 2 O 
x (cm) x (cm) 
Figure 11. Axisymmetric compression corner. Surface profiles. 
(a) Static pressure. (b) Skin friction coefficient. 
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Table 3. Axisymmetric compression corner. Separation and reattachment locations. 


Separation Reattachment Separation 

Location (cm) Location (cm) Length (cm) 
2: snnanegneee -2.730 0.970 3.700 
(approximate values) 


The axial velocity profiles are shown in Figure 12. At the first measurement station, at x = -4.5 cm, the flow is 
attached and the velocity profile thickness of the computations matches the experiment, with the exception of the 
SSG/LRR result, which is separated. The SA, SST-V, and Wilcox models agreed well with each other and with the 
data, though they were slightly less full at the wall. At x = -3.0 cm, the solutions computed with the SST-V and 
Wilcox models are separated, and the SA model solution and the experimental boundary layers remain attached. All 
of the models show separation at x = -2.0 cm, where the experiment is also separated. All of the models give a stronger 
and larger separated region than the experiment. The SA model result is the closest to the experimental value; the 
SST-V and Wilcox models behave very similarly and show a stronger and larger separation, and the SSG/LRR RSM 
drastically over predicts the strength of the separation resulting in a very different profile shape. Throughout the 
separated region the SST-V and Wilcox models behave similarly. In the reattached region on the flair, at x = 3.464 
cm, the two models have different profiles. However, once the flow attaches, all of the profiles approach the 
experiment, as shown at x = 5.896 cm, with the SSG/LRR result still deviating the most. 
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Figure 12. Axisymmetric compression corner. x-velocity profiles. 


The turbulent shear stress profiles are shown in Figure 13. The SA and SSTV results are similar to those computed 
with Wind-US and reported by DeBonis [12]. At the most upstream station, x = -4.5 cm, the SA, SST-V and WRSM 
models under-predict the shear stress near the wall, but agree with the experiment in the outer portion of the boundary 
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layer. The SSG/LRR RSM solution is separated at this station, and the profile somewhat resembles the experimental 


profiles further downstream in the separated region at x = -2.0 cm. Overall, the CFD shear stress profiles are in poor 
agreement with the experiment at all stations. 


a5 x =-4,500 cm =-3.000 cm 
Experiment 
—— SA 
| ——— Sst 
———— SSG/LRR 
1-}|————._ Wilcox RSM 


R004 "0.002 ; 0.002 8.07 "0.005 , 0 0.005 
<uv>/uy <uv>/u) 
= -2.000 cm x =-1.000 cm 
25 
2 
15 
£ 
2 
P 
>, 7 
6 
5 
: o 
0.5 = 5 
3 
8a 9.005 =O 0.005 8.01 "0.005 , (0 0.005 
<uv>/u> <uv>/u>, 
x = 0.433 cm x = 1.732 cm 


2 
Poi 5 0.01 0.005 0 0.005 Poi 5 0.01 0.005 , 0 0.005 
<uv>/u" <uv>/u) 
x = 3.464 cm x = 5.896 cm 
3.5 3.5 
3 3 
2.5 25 
e? g° 
2 2& 
>1.5 >1.5 


° 
R006 0.004 = -0.002 ? 0.002 0.004 
<uv>/u_, 


-f006 -0.004 0.002 2 0.002 0.004 
<uv>/u 


Figure 13. Axisymmetric compression corner. Turbulent shear stress profiles. 
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Overall, for the supersonic axisymmetric compression corner, the SA model does the best job of computing the 
onset of separation and the velocity profiles. The Wilcox RSM is very similar to the SST-V model, and the SSG/LRR 
RSM does a poor job of computing this flow. For this case, use of the RSMs was found to provide no benefit. 


C. Mixing Layer 

1, Experimental Configuration 

In this case, two supersonic streams, initially separated by a splitter plate, come into contact and form a shear layer. 
Results are compared with test case 2, in the experiment of Goebel and Dutton [14]. At the entrance of the mixing 
layer, the flow conditions are M,=1.91, T,,=578 K, and U,;=700 m/s for the high-speed stream and M,=1.36, T;,,, =295 
K, and U,=399 m/s for the low-speed stream. Both streams have a pressure of 49 kPa. The relative Mach number, 
M, = AU/(@) =0.91, where AU = (U, — U,), and @ is the average of the freestream speed of sound. 


2. CFD Methodology 

The Wind-US computations shown below are nearly identical to those reported by Yoder [20], with only some 
minor differences in the grid. The 2-dimensional structured grid consisted of three zones. Zones 1 and 2 were 
upstream of the mixing layer, each having 101 points in the axial direction and 89 points in the vertical direction. The 
lengths of the incoming streams were chosen to match the experimental values of the momentum thickness. Lengths 
of 240 mm and 160 mm were required for the high-speed and low-speed streams, respectively. The third zone, 
beginning at the end of the splitter plate was 401x241 and 500 mm in the axial direction by 48 mm in the vertical 
direction. The plate was 0.5 mm thick and had 64 points along the vertical edge. It is shown in Figure 14. The grid 
was clustered vertically on the splitter plate to capture the boundary layers with wall spacing of 0.0025 mm resulting 
in y* values of 0.8-1.0. It was also clustered axially near the splitter plate trailing edge such that x* was approximately 
16. Results computed using a version of this grid with every other point removed were nearly equivalent to results 
computed using the full grid. Viscous adiabatic wall boundary conditions were used on the splitter plate. In the 
experiment, the top and bottom test section walls were angled to account for boundary layer growth; in the 
computations, inviscid boundary conditions were used. At the inflow, the flow was fixed, and at the downstream 
boundary, the flow was extrapolated. 

The grid used for FUN3D was created from the Wind-US grid, with a few modifications. It was extruded 1mm in 
the z-direction, to make it 3-dimensional. The tight axial packing at the inflow boundary was not compatible with 
FUN3D, so a short region (10mm) with inviscid walls was added upstream of the viscous walls and the axial clustering 
was relaxed. The total number of axial grid points remained the same. 

For the Wind-US computations, the SST-V model and an Explicit Algebraic Stress Model (EASM) [20] were 
used. The EASM is a k-epsilon formulation that uses a nonlinear equation to compute the shear stress tensor, and is 
described in Ref. [20]. Results computed using this model are included for comparison because it gives a more 
complete calculation of the shear stress tensor compared to models based on the Boussinesq approximation while 
remaining less computationally expensive than the RSMs. For the FUN3D computations, the SST-V model, the 
SSG/LRR, and Wilcox RSMs were used. 


it) 
“n 
-200 -100 it) 100 x{mm) 200 300 400 600 
Figure 14. Mixing layer. Computational mesh. 
(Every fourth grid point shown.) 


3. Results 

Mean velocity, turbulent intensity, and turbulent shear stress profiles at several axial stations are given in Figure 
15-Figure 18, and the shear layer thickness is given in Figure 19. The Wind-US SST-V results are also included in 
these plots and are in excellent agreement with the FUN3D SST-V results. At the entrance to the mixing section, 
Figure 15-Figure 18 show that all of the turbulence models give nearly the same profile and agree well with data, 
however downstream, the results deviate. The velocity profiles are shown in Figure 15 and indicate that, beginning 
at x = 100 mm, the SST-V model has a linear slope in the mixing region, with a sharp change to the freestream values 
at both edges. The EASM model has a nonlinear transition at the edges of the mixing region, and the Wilcox and 
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SSG/LRR RSMs behave nearly the same and their profiles fall somewhat in between the SST-V and EASM profiles 
for the curvature at the mixing-layer edges. The streamwise turbulence intensity is shown in Figure 16. At all axial 
stations, the SST-V model under-predicts the maximum values, while the EASM does well at computing the maximum 
values in the center of the mixing section. The SSG/LRR model predicts similar maximum values, while the Wilcox 
RSM over-predicts the maximum values. The profiles are fairly similar in shape, with the EASM model having the 
most gradual transition at the edges of the mixing region. The transverse turbulence intensity is given in Figure 17. 
Again, at the entrance to the mixing section, the models all give essentially the same result. Moving downstream, the 
SST-V model significantly over-predicts the maximum values. The EASM and the RSMs do a good job predicting 
the peak values downstream, with the EASM having a more gradual slope at the edges of the mixing region and giving 
a slightly more narrow mixing region. The turbulent shear stress is shown in Figure 18. The EASM does the best job 
of predicting the maximum values and the shape. The two RSMs give very similar results with slightly higher 
maximum values; the SST-V results are very similar, with a slightly wider mixing region. 

The shear layer thickness is given in Figure 19 and is defined as the distance b, between transverse locations where 
U =U, — 0.1AU and U = U,+0.1AU. The RSMs do the best job at getting the shear layer thickness correct. The 
SST-V model results have the highest values and the EASM gets the downstream values too low. This corresponds 
to the streamwise velocity and turbulence quantity profiles in Figure 15-Figure 18. The width of the shear layer 
predicted by the SST-V is largest, with the RSMs slightly smaller. The RSMs have slightly more curvature at the 
edges of the shear layer than the SST-V model and less than the EASM. Since the definition of b uses the 10%AU 
criteria, the curvature at the edges of the shear layer effects the value. 

While all of the turbulence models examined appear to predict the streamwise velocity profiles well, it appears 
that the EASM and the RSMs have benefits over the two-equation SST-V model. Overall, they are closer to the 
experimental values in predicting the turbulence intensity, turbulent shear stress and shear layer thickness. This is 
most likely due to the fact that they account for anisotropy effects in their calculation of the shear stress tensor. 
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Figure 15. Mixing layer. Mean velocity profiles. 
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Figure 16. Mixing layer. Streamwise turbulence intensity. 
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Figure 17. Mixing layer. Transverse turbulence intensity. 
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Figure 18. Mixing layer. Turbulent shear stress profiles. 
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Figure 19. Mixing layer. Shear layer thickness. 


D. Axisymmetric Jet 

1, Experimental Configuration 

In this experiment, subsonic, unheated air is passed through an axisymmetric convergent nozzle to produce a Mach 
number of Mjer = Ujer/Ajer = 0.51, exiting into quiescent air. (The acoustic Mach number is, Mg = Ujer/Grer = 
0.5.) This experiment used an acoustic research nozzle referred to as ARN2, having an exit diameter of 2 inches, at 
conditions referred to as set point 3. [35] The freestream static temperature and pressure were 530 degrees R and 14.3 
psi, respectively. Particle image velocimetry (PIV) data was taken of axial and turbulent velocities. 

2. CFD Methodology 
This case was computed using the Wind-US code with the SA and SST-V turbulence models as described on the TMR 
[19]. Since quiescent air is often difficult to compute with CFD, the flow external to the nozzle was set to Mer = 
0.01, where the external flow conditions are referred to as the reference conditions. At the jet inflow, 
Pt/Pre¢ = 1.197 and T,/T,¢¢ = 1.0. The nozzle walls had adiabatic viscous boundary conditions and the jet outflow 
was extrapolated. A series of grids was available on the TMR, each coarse grid was created by removing every other 
grid point from the next finest grid. The grids contained 3 zones each as shown in Figure 20 and Table 4. To convert 
these grids to FUN3D-compatible format, Pointwise was used, resulting in a single zone grid. The grids were one- 
cell deep spanning approximately | degree in the theta direction. For computations using the RSMs, a 90-degree sector 
was required. To generate this grid, the axisymmetric grid was extruded 90 degrees in the theta direction at 5-degree 
increments for a total of 19 grid points. A grid resolution study was done using FUN3D with the SA and SST-V 
turbulence models. The solution did not show appreciable changes when using Grids 129 and Grid 257, so grid 129 
was used for the remainder of the computations. 

Regardless of the turbulence model chosen, none of the solutions converged to a definite steady state result. When 
run in a time-accurate mode, solutions using the SA, SST-V and SSG/LRR model converged to a quasi-steady result 
when the solutions were time-averaged, however the solution computed using the Wilcox RSM was still noticeably 
changing after over one million time steps, and so the results are not included in the plots. The solutions presented on 
the TMR were computed with the CFL3D [36] and the Wind-US codes, and also required a time-accurate solution 
method to achieve a quasi-steady state result. Wind-US results using the SA and SST-V turbulence models on Grid 
257, were also given on the TMR. FUN3D results computed on Grid 257 computed using the SA and SST-V models 
were in good agreement with these Wind-US results. 
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Table 4. Axisymmetric jet. Grid study. 


Name Upstream External Region Internal Nozzle Region Jet Mixing Region 
Grid 65 29529 16x25 65x57 
Grid 129 49x49 31x49 129x113 
Grid 257 07x97 61x97 257x225 
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Figure 20. Axisymmetric jet. Grid shown at the second from the finest level (grid 129). 


3. Results 

The results are given in Figure 21-Figure 26 below and include centerline velocity and turbulent kinetic energy 
profiles, and radial velocity, turbulent shear stress and turbulent kinetic energy profiles at 5 locations downstream of 
the jet exit. As mentioned above, the solution computed using the Wilcox RSM did not achieve a quasi-steady result, 
and the results are omitted from the plots. The mean x-velocity is given in Figure 21 and indicates that the SSG/LRR 
RSM does the best job of predicting the overall mixing. Although the onset of mixing takes effect at about 6.5 
diameters downstream of the jet exit, 2 nozzle-exit diameters downstream of experimental onset location of x/Djer = 
4.5, the SSG/LRR RSM result otherwise gives the best overall mixing throughout the mixing region. The SA model 
result shows the onset of mixing closest to the experiment at x/Dj., = 5.8, however it over-mixes, resulting in a lower 
downstream centerline velocity. The SST-V model result shows the onset beginning the furthest downstream of the 
experiment at x/Dje, = 8.1, then over-mixes, resulting in downstream velocities similar to the SA result. The 
turbulent kinetic energy, where k for the SA result was computed using the Boussinesq approximation and Bradshaw’s 
approximation as described in Section IV.A.3, is plotted in Figure 22. The general profile shape is similar to that of 
the experiment. However, the values begin to increase further downstream, and the peak values are quite different. 
The SA model peak turbulent kinetic energy is 60% lower than the experiment, the SST-V peak value is 22% higher 
and the SSG/LRR peak is 42% higher. All of the CFD results show that the turbulent kinetic energy begins to increase 
between 5 and 8 jet diameters downstream of the jet exit, whereas the experimental turbulent kinetic energy begins to 
rise at the jet exit. For the SA model, the slope as the turbulent kinetic energy increases is closest to the experiment; 
the SST-V and SSG/LRR RSM have much steeper slopes. 
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Figure 21. Axisymmetric jet. Axial velocity profiles on the centerline computed using FUN3D. 
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Figure 22. Axisymmetric jet. Turbulent kinetic energy profiles along the centerline computed using FUN3D. 


The x-velocity profiles are shown in Figure 23 at x/Djet = 2.,5.,10.,15.,and 20. Near the jet exit, all of the 
profiles agree well with the experiment and predict the correct profile shapes. Beginning at x/Dje¢ = 10., the values 
near the centerline are less than the experiment for the SA and SST-V computations, as was shown in Figure 21, but 
agree well away from the centerline. The y-velocity, turbulent shear stress and turbulent kinetic energy are shown in 
Figure 24-Figure 26, and show how the computations predict the shear layer mixing. At x/Dje¢ = 2. and 5., the 
SSG/LRR RSM does slightly better at predicting the profile shape. However, further downstream, the y-velocity and 
the turbulent shear stress results are similar for all three models. At the three most downstream stations, the turbulent 
kinetic energy is different at the centerline, as was also indicated by the centerline profile in Figure 22. Overall for 


this subsonic axisymmetric jet study, the SSG/LRR model shows some benefits over the SA and SST-V models at 
predicting the mixing. 


25 
American Institute of Aeronautics and Astronautics 


Experiment 
SSG/LRR RSM 


SA 
SST 


te) 
N 
=) 


Figure 23. Axisymmetric jet. Radial x-velocity profiles computed using FUN3D. (Subsequent profiles shifted by 
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Figure 24. Axisymmetric jet. Radial y-velocity profiles computed using FUN3D. (Subsequent profiles shifted 
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Figure 25. Axisymmetric jet. Radial turbulent shear stress profiles computed using FUN3D. (Subsequent 
profiles shifted by u'v'/Ujer” = 0.01.) 
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Figure 26. Axisymmetric jet. Radial turbulent kinetic energy profiles computed using FUN3D. (Subsequent 
profiles shifted by k/Ujer” = 0.03.) 


E. Conclusion 

Two full second-moment Reynolds stress turbulence models available in the FUN3D code, the Wilcox and the 
SSG/LRR, were evaluated for four test cases: a transonic diffuser, a supersonic axisymmetric compression corner, a 
supersonic compressible planar mixing layer, and a subsonic axisymmetric jet. These model results were compared 
with solutions computed using the SA and SST-V turbulence models, and for the planar mixing layer case, an EASM. 
The transonic diffuser was the first case evaluated, and gave inconclusive results as to the benefits of the model. The 
supersonic axisymmetric compression corner, which has a complex shock system and a large region of separated flow, 
was evaluated next. The SA model did the best at computing the pressure rise and the separation location and length, 
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and the Wilcox RSM gave results similar to the SST-V model. The SSG/LRR RSM solution had oscillatory residuals, 
and severely over-predicted the onset of separation. All of the models had difficulty computing the boundary layer 
profiles and turbulence quantities in the separated region, and no additional benefit was provided by using the RSMs 
for this case. For the supersonic planar mixing layer, both RSMs and the EASM had benefits over the SST-V model 
at predicting the turbulence intensity, turbulent shear stress and shear layer thickness. These benefits were most likely 
because they account for anisotropy effects in their calculation of the shear stress tensor. For the subsonic 
axisymmetric jet, the SSG/LRR model predicted the mixing of the core velocity the best, whereas the solution with 
the Wilcox RSM did not reach a stable answer for this case. 

The four cases examined are flows that are challenging for current turbulence models because they contain mixing, 
shock waves and/or separation. Overall, the RSMs showed benefit over the SA and SST-V models for the planar 
mixing layer and the axisymmetric jet flow, and may be useful for future nozzle calculations. The EASM, which is 
less cumbersome to code and requires less computational time, did well with the planar mixing layer, and also may 
be a viable choice for mixing flows. While the cases examined are challenging flows, they are still relatively simple 
in geometry and flow features. More complex flow cases may reveal more benefits of the RSMs and are recommended 
for future study. 
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